version 19
drop _all

* change the root folder according to your computer directory
global root = "D:\WorkingPaper-Series\Township_Light\Replication_Files"  

global dofiles = "$root\dofile"     
global working_data = "$root\workdata"
global tables = "$root\table"
global figures = "$root\figure"

cd "$working_data"


******************
*****Pre Set******
******************

global geo_ctrl = "topography_t* slope_t* roughness_t*"
global wth_ctrl = "gaez_t* precipitation_t* sunlight_t*"
global dist_ctrl = "county_dist_t* dist_hz_t* dist_coast_t*"

global full_ctrl = "$geo_ctrl $wth_ctrl $dist_ctrl"

global geo_ctrl_1 = "topography_t_1 slope_t_1 roughness_t_1"
global wth_ctrl_1 = "gaez_t_1 precipitation_t_1 sunlight_t_1"
global dist_ctrl_1 = "county_dist_t_1 dist_hz_t_1 dist_coast_t_1"

global full_ctrl_1 = "$geo_ctrl_1 $wth_ctrl_1 $dist_ctrl_1"

global geo_ctrl_2 = "topography_t_2 slope_t_2 roughness_t_2"
global wth_ctrl_2 = "gaez_t_2 precipitation_t_2 sunlight_t_2"
global dist_ctrl_2 = "county_dist_t_2 dist_hz_t_2 dist_coast_t_2"

global full_ctrl_2 = "$geo_ctrl_2 $wth_ctrl_2 $dist_ctrl_2"

global matchvar = "county_dist topography slope gaez sunlight dist_hz dist_coast"


global pre_ctrl = "light_t road_t pop_t"


*********************
***** Figure C5 *****
*********************


use "township_panel_final.dta", clear

collapse (mean) countycode citycode first_year treat light10 road13 population10 county_dist topography digit_village rural_revital total_freq, by(idcode)

gen log_pop = log(population10)
gen log_topo = log(topography+1)
gen log_road = log(road13)
replace county_dist = county_dist/10
preserve
***** Before and After 2016
foreach w in county_dist light10 log_pop log_road log_topo {
	mat define est_`w' = J(2,2,0)
	qui reghdfe `w' treat if first_year >=2016|first_year ==0, a(citycode) cl(citycode)
	mat est_`w'[1,1] = e(b)[1,1]
	mat est_`w'[1,2] = sqrt(e(V)[1,1])
	qui reghdfe `w' treat if first_year <2016|first_year ==0, a(citycode) cl(citycode)
	mat est_`w'[2,1] = e(b)[1,1]
	mat est_`w'[2,2] = sqrt(e(V)[1,1])
}

foreach w in county_dist light10 log_pop log_road log_topo {
	svmat est_`w'
	rename (est_`w'1 est_`w'2) (coef`w' se`w')
}

keep coef* se*
keep in 1/2
gen after_2016 = -1*_n+2

reshape long coef se, i(after_2016) j(type) string

export delimited using "feature_after_2016.csv", replace
restore

preserve
***** High KW versus low KW
foreach w in county_dist light10 log_pop log_road log_topo {
	mat define est_`w' = J(2,2,0)
	qui reghdfe `w' treat if total_freq>=5, a(citycode) cl(citycode)
	mat est_`w'[1,1] = e(b)[1,1]
	mat est_`w'[1,2] = sqrt(e(V)[1,1])
	qui reghdfe `w' treat if total_freq<5, a(citycode) cl(citycode)
	mat est_`w'[2,1] = e(b)[1,1]
	mat est_`w'[2,2] = sqrt(e(V)[1,1])
}

foreach w in county_dist light10 log_pop log_road log_topo {
	svmat est_`w'
	rename (est_`w'1 est_`w'2) (coef`w' se`w')
}

keep coef* se*
keep in 1/2
gen high_kw = -1*_n+2

reshape long coef se, i(high_kw) j(type) string

export delimited using "feature_high_kw.csv", replace
restore

